function [ownElas, derMat] = rcnl_der( ...
    alpha,         ...
    Pi,            ...
    rho,           ...
    s_jt,          ...
    p_jt,          ...
    xnonlin,       ...
    cdindex,       ...
    dfull,         ...
    cdid,          ...
    prodid,        ...
    delta,         ...
    expmval_mat,   ...
    replace,       ...
    rho_transform, ...
    tolerance      ...
)

if rho_transform
    rho = rho_transform * (1/(1+exp(-rho)));
end

if isnan(delta)
    [delta, expmu] = rcnl_meanval(Pi, rho, s_jt, xnonlin, cdindex, dfull, expmval_mat, replace, tolerance);
else
    expmu = ones(size(dfull));
end
pcoefi = alpha + dfull .* Pi(1);

[sharei, scondi, ~] = rcnl_indsh(exp(delta), expmu, rho, cdid);
numProd = length(unique(prodid));
derMat  = zeros(numProd, numProd, max(cdid));
ownElas = zeros(size(s_jt));

for m = 1:max(cdid)
    sel = (cdid == m);

    shareiM  = sharei(sel, :);
    scondiM  = scondi(sel, :);
    pcoefiM  = pcoefi(sel, :);

    der     = rcnl_der_price(rho, shareiM, scondiM, pcoefiM);
    emat    = der .* repmat(p_jt(cdid==m)', [size(der, 1) 1]) ./ repmat(s_jt(cdid==m), [1 size(der, 1)]);

    derMat(prodid(sel), prodid(sel), m) = der;
    ownElas(sel) = diag(emat);
    
end
